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1     INTRODUCTION 

All  objects  in  an  atmospheric  environment  have  a  potential  associated  with  them    This 
charge  is  usually  much  larger  for  moving  airborne  objects  than  for  stationary  (ground  contacting) 
objects.  This  report  presents  preliminary  results  obtained  during  Fall  1993  in  the  area  of 
Electrostatic  Target  Detection.  This  preliminary  study  investigates  the  application  of  time-domain 
filtering  methods  to  estimate  the  "optimum"  fuzing  location  and  studies  the  robustness  of  the 
method  to  white  noise  degradation 

Section  2  presents  the  problem  set-up  Section  3  describes  the  methodology  investigated 
during  this  preliminary  study,  its  robustness  in  the  presence  of  additive  white  noise  and  presents  a 
methodology  for  colored  noise  generation.  Finally,  Section  4  presents  recommendations  for 
further  research. 
2.  PROBLEM  SET-UP 

Powered  and  moving  objects  tend  to  charge  electrostatically  due  to  combination  of 
different  charging  processes.  This  charge  is  usually  much  larger  for  airborne  objects  than  for 
ground  contacting  objects.  Measurements  and  some  preliminary  experiments  have  shown  that, 
under  ideal  conditions,  one  can  sense  the  presence  of  objects  of  interest  using  electrostatic  sensors 
(ESS).   In  this  preliminary  study  we  use  the  point  charge  response  of  a  short-circuit  longitudinal 
sensor  to  enhance  the  detectability  of  moving  objects,  that  is  detection  prior  to  passing  the  object 
under  ideal  and  noisy  conditions  [1].  The  information  available  for  the  proximity  fusing  sensing 
scheme  is  the  differentiated  field  in  the  direction  of  the  projectile.   An  example  for  a  typical 
Electrostatic  point-charge  response  is  given  in  Figure  1,  and  further  details  regarding  its 
expression  can  be  found  in  Section  3.   One  key  feature  was  determined  to  be  the  position  of  the 
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zero-crossing  [1].   As  a  result,  we  concentrated  in  estimating  the  location  of  the  zero-crossing  in 
this  preliminary  study,  and  investigated  the  robustness  of  the  technique  to  noise  degradation 

The  main  motivation  behind  the  choice  of  the  technique  investigated  in  the  preliminary 
study  lies  in  its  simplicity,  as  it  can  be  implemented  real-time  and  with  a  minimum  of  hardware 
Estimating  the  position  of  the  zero-crossing  when  no  noise  is  present  is  a  very  basic  problem 
However,  for  practical  applications  one  has  to  deal  with  noisy  responses,  such  as  given  in 
Figure  2    In  such  a  case,  the  problem  becomes  to  estimate  the  location  of  the  "correct" 
zero-crossing  located  around  x  =  -5.653  1  meters.   Note  that  due  to  noise  degradations,  numerous 
zero-crossing  points  can  be  observed  before  the  optimum  fuzing  point    As  a  result,  we  use  an 
averaging  filter  to  smooth  the  effects  of  the  noise,  and  to  insure  that  there  is  no  zero-crossing 
before  the  optimum  fuzing  point.   However,  as  a  consequence,  the  averaging  filter  flattens  out  the 
response  and  introduces  further  error  in  the  location  of  the  estimated  fuzing  point    This  problem 
can  be  minimized  by  using  a  non-linear  median  smoothing  filter.  The  median  filter  preserves  sharp 
and  sustained  changes  in  the  response,  and  as  such  preserves  the  slope  of  the  response  better  than 
the  averaging  filter  does  while  smoothing  some  of  the  noise-like  components  out    However,  the 
median  filter  usually  does  not  smooth  the  noise-like  components  sufficiently    As  a  result,  a  good 
solution  may  be  obtained  with  an  algorithm  which  has  a  combination  of  first  an  averaging  filter, 
and  a  second  median  filter,  to  insure  a  better  estimation  of  the  fuzing  point    Mean  and  median 
filters,  and  the  definition  of  the  threshold  needed  to  switch  from  one  filter  type  to  the  other  are 
discussed  in  further  details  in  Section  3. 


3    METHODOLOGY  INVESTIGATED 

3.1.   INTRODUCTION  TO  THE  MEAN  AND  THE  MEDIAN  FILTER 
3.1.  a.  MEAN  FILTER 

The  mean  filter  is  an  averaging  process  that  acts  as  a  smoothing  filter  as  the 
technique  outputs  the  average  of  the  input  points  considered.  Note  that  only  causal  outputs  are 
used  for  the  purpose  of  this  study.   A  causal  system  is  a  nonpredictive  system  which  means  the 
output  does  not  precede  the  input  [2]    The  causal  mean  function  (f(al)  expressed  as  a  function  of 
the  input  function  (f(j))  is  given  by: 

f(a)  (n  +  L  - 1 )  =  mean(f(l)(n  to  n  +  L  -  1 )), 
where  n  is  the  index  of  the  data  points  and  L  is  the  length  of  the  input  sequence  being  considered 
by  the  filter.   Figure  3  and  4  show  an  input  and  the  result  of  smoothing  the  input  by  using  a  mean 
filter  of  length  three  (L  =  3). 

3.1.  b.  MEDIAN  FILTER 

The  median  filter  is  a  non  linear  filter  that  outputs  the  median  of  the  input 
points  considered.  Note  that  as  a  result  of  the  non  linearity  structure,  the  slope  of  the  input  curve 
is  preserved  better  than  when  the  mean  filter  is  used  [3]    The  median  filter  is  used  in  a  causal 
configuration,  and  the  median  function  f(m)  is  expressed  as  a  function  of  the  input  function  f(i, 
given  by: 

f(m)(n  +  L  -  1)  =  median(f(u(n  to  n  +  L  -  1)). 
Figure  3  and  4  show  an  input  and  the  result  of  smoothing  the  input  by  using  a  median  filter  of 
length  three. 


3.2.  MODEL 


The  model  used  here  is  the  point-charge  response  of  a  short-circuit  longitudinal 
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sensor,  where  the  electrostatic  point-charge  response  '—f-  has  been  derived  [1]  to  be 


dEn  _    a  (     y~-2x- 


dt 


dt         V(.t:+r)5/:: 
where  A  =  j^-  and  Q  is  the  charge  (in  Coulombs),  vxis  the  probe  velocity  (in  the  x-direction,  in 

m/s),  x  is  the  distance  between  the  probe  and  the  point  of  closest  approach,  and  y  is  the  miss 
distance  (orthogonal  to  x,  in  meters).  Figure  1  is  an  illustration  of  the  normalized  point  charge 
response  with  A  =  1  and  y  =  8  meters 

The  portion  of  the  point  charge  response  we  concentrate  on  is  the  part  before  the 
explosion  point  is  reached,  or  in  this  case,  the  part  before  the  first  zero  crossing,  as  shown  in 
Figure  5.   The  curve  is  composed  only  of  very  small  frequency  components,  as  shown  by  its 
Fourier  transform  in  Figure  6.   As  a  result,  any  low  pass  filter  cancels  the  effects  of  a  colored 
noise  centered  outside  of  the  model  bandwidth,  illustrated  in  Figures  7  and  7. 1 .  Thus,  the  most 
critical  type  of  noise  to  consider  is  the  white  noise  case  as  it  covers  the  entire  frequency  spectrum 
Figure  2  shows  the  point  charge  response  with  additive  white  noise  with  a  noise  variance  equal  to 
1/49. 

3  3    DETECTION  TECHNIQUE 

3  3a    DETECTION  OF  THE  FUZING  POINT  -  NOISE-FREE  CASE 

The  fuzing  point  is  located  at  the  first  zero-crossing  of  the  point  charge 
response,  as  illustrated  in  Figure  5.  This  first  zero-crossing  is  computed  by  tracking  the  change  of 
sign  of  the  point  charge  response.  The  fuzing  point  value  is  found  to  be  equal  to  -5.653  1  meters 
in  the  x  direction. 


3.3  b    DETECTION  OF  THE  FUZING  POINT  -  NOISY  CASE 

When  additive  noise  is  present,  the  above  simple  detection  procedure  is  no 

longer  valid,  as  numerous  zero-crossings  occur  due  to  the  noise  before  the  correct  fuzing  point 

value.  The  effects  of  the  noise  can  be  minimized  by  applying  a  simple  low-pass  averaging  filter. 

This  filter  has  the  effect  of  suppressing  early  zero-crossings.  However,  it  also  has  the  drawback 

of  attenuating  the  slope  in  the  region  where  the  true  fuzing  point  occurs,  and  introduces  further 

error  in  the  estimation  of  the  fuzing  point    A  non-linear  filter,  such  as  the  median  filter  introduced 

in  Section  3.1.b,  can  be  used  to  preserve  sudden  sustained  value  changes  in  the  data  and  filter 

some  noise  components    However,  a  median  filter  does  not  provide  enough  smoothing  of 

noise-like  components  where  the  noisy  point  charge  response  is  very  small    As  a  result,  a  good 

compromise  may  be  reached  by  using  a  combination  of  both  mean  and  median  filters  to  minimize 

the  effects  of  the  noise  and  better  estimate  the  location  of  the  fuzing  point.  The  algorithm  starts 

with  a  mean  filter  to  smooth  the  effects  of  the  noise  down  and  to  insure  that  no  early 

zero-crossings  occur.   Once  the  filtered  point  charge  response  values  fall  below  a  predefined 

threshold,  set  by  the  user,  the  filtered  point-charge  response  is  insured  not  to  cross  the  zero  axis 

early    When  the  threshold  value  is  reached,  the  detection  procedure  switches  to  the  median  filter 

The  median  filter  still  smoothes  the  effects  of  the  noise  down  and  better  preserves  the  slope  of  the 

point-charge  response  in  the  fuzing  point  neighborhood.  Estimation  of  the  fuzing  point  location  is 

determined  when  the  noisy  response  changes  sign    Values  for  the  length  of  the  mean  filter,  the 

median  filter,  and  the  threshold  value,  are  all  set  by  the  user.  The  resulting  point  charge  response 

is  defined  as: 

s„(k)  =  s(k)  +  oww(k), 


where  s(k)  is  the  noise-free  normalized  point  charge  response,  w(k)  is  an  additive  normal  white 
noise  with  mean  equal  to  zero  and  variance  equal  to  one    The  parameter  aM  is  adjusted  during  the 
calculations  to  vary  the  power  of  the  noise    The  resulting  Signal-to-Noise  Ratio  (SNR)  for  the 

noisy  point  charge  response  is  defined  in  decibels  as 

-101og(4), 

where  s  is  the  point  charge  response  power  between  -40  and  40  meters.  Figures  8  through  9. 1 

show  the  noisy  point  charge  response  and  the  filtered  noisy  point  charge  response  with  the  change 
of  filters  indicated  by  an  asterisk  (*)  for  three  additive  noise  power  levels 

The  threshold  equal  to  -0.20  and  length  of  the  filters  equal  to  15  are 
determined  to  insure  a  reliability  of  99%  .   The  reliability  measure  [4]  is  calculated  in  the 


following  fashion: 


R=] 


N    ' 

where  R  is  the  reliability,  n  is  the  number  of  failures,  and  N  is  the  number  of  trials    Figure  10 

shows  the  normalized  error  measure  defined  as: 

tr  =  — ; , 

where  eesl  is  the  estimated  zero  crossing  location  obtained  using  the  proposed  technique,  and 

elh  =  -5.653  1  is  the  noise-free  zero  crossing  location  in  meters,  for  340  trials  of  the  combination 

mean/median  filtered  noisy  point  charge  response    In  Figure  10  the  noise  variance  is  equal  to 

1/100,  i.e  the  SNR  is  equal  to  8.2086  dB,  and  there  is  no  error  greater  than  one  meter  from  the 

fuzing  point  which  yields  a  reliability  greater  than  99%    The  parameters  for  the  filters  are  length 

(L)  and  the  threshold  (E)  set  to  L  =  1 5  and  E  =  -0.2. 

3  4    ROBUSTNESS  OF  MEAN/MEDIAN  FILTER  COMBINATION 

The  robustness  of  the  mean/median  filter  combination  is  studied  by  varying  the 

power  of  the  white  noise  added  to  the  model    One  hundred  ( 1 00)  trials  are  computed  for  each 


white  noise  power.  Figure  1 1  shows  the  mean  and  two  standard  deviations  of  the  normalized 
error    The  results  show  that  the  mean  normalized  error  is  maintained  to  within  a  meter  for  a  noise 
variance  less  than  1/55,  SNR  less  than  5.6122  dB.  The  standard  deviations  for  two  SNR 
realizations  are  near  0.2,  which  is  the  result  of  a  trial  at  that  noise  variance  having  an  early  zero 
crossing    Further  improvement  for  these  two  SNR  realizations  can  be  obtained  by  decreasing  the 
value  of  the  threshold  to  avoid  an  early  zero  crossing. 

3.5.  COLORED  NOISE  MODELS 

Colored  noise,  formed  as  band  limited  Gaussian  noise,  can  be  used  to  simulate 
natural  known  noises  in  the  target  area.   The  MATLAB  code  implementation  in  the  Appendix  is 
menu  driven  and  can  be  set  to  prompt  the  user  for  center  frequency  (fc)  and  a  quality  factor  (Q): 

Q  =  fc/(fu  -  fl), 
where  fu  is  the  upper  cutoff  frequency  and  fl  is  the  lower  cutoff  frequency.  The  program  can  also 
be  set  to  accept  the  cutoff  frequencies  directly  along  with  the  sampling  frequency.   The  sampling 
frequency  is  important  to  place  the  band  of  the  noise  in  the  proper  place  on  the  frequency 
spectrum    Figures  12  through  14.1  present  examples  of  colored  noise  realizations  and  their 
spectra  with  different  Q  values. 
4    RECOMMENDATIONS  FOR  FURTHER  RESEARCH 

This  report  investigated  the  applicability  of  time-domain  techniques  to  estimate  the  fuzing 
location  in  the  electrostatic  detection  problem    This  preliminary  study  is  restricted  to  the  case  of 
an  ideal  point  charge  response  for  a  short  circuit  longitudinal  sensor,  and  a  miss  distance  fixed  to 
eight  meters.  The  main  advantage  of  this  analysis  technique  is  its  simplicity;  a  minimum  of 
software  would  be  needed  for  its  implementation    The  procedure  uses  a  combination  of  a  linear 


smoothing  filter  and  a  non-linear  median  filter  to  estimate  the  fuzing  location  in  noisy  conditions 
The  robustness  of  the  technique  to  additive  white  noise  distortions  is  also  investigated    Initial 
results  show  that  the  time-domain  analysis  technique  is  promising;  the  error  obtained  in  the 
estimation  my  be  kept  to  within  10%  of  the  true  fuzing  location  for  a  Signal-to-Noise  Ratio 
(SNR)  greater  than  5.4  IdB  (i.e.,  for  an  additive  noise  power  less  than  0.02)  in  86%  of  the 
distortion  levels  investigated.  For  two  of  the  SNR  realizations  investigated  the  error  goes  up  to  a 
maximum  of  27%  of  the  true  value.  Further  improvements  can  be  obtained  by  refining  the 
choices  for  the  threshold,  and  the  filter  lengths 

Several  points,  which  have  not  been  considered  in  this  preliminary  study,  need  to  be 
addressed  if  the  study  is  to  be  taken  any  further    One  needs  to  investigate  the  robustness  of  the 
procedure  to  specific  additive  colored  noise  distortions.  Note  that  due  to  the  low  frequency 
contents  of  the  ideal  point-charge  response,  only  colored  noise  which  have  bandwidths  that 
overlap  the  bandwidths  of  the  point  charge  response  need  to  be  investigated    Furthermore,  we 
have  restricted  our  study  to  a  time-domain  analysis.   Improvements  may  potentially  be  obtained 
by  investigating  frequency  techniques  and  correlation  procedures,  and  further  work  is  needed  to 
compare  their  performances.   In  addition,  the  results  presented  in  this  preliminary  study  are 
restricted  to  the  case  of  a  point  charge  response  for  a  longitudinal  sensor,  and  for  a  fixed  miss 
distance.  Typical  signatures  are  more  complex  [1],  and  one  needs  to  assess  the  applicability  of 
any  technique  investigated  to  real-time  signatures 


Figure  1   Normalized  point  charge  response  with  the  miss  distance  equal  to  eight 
meters. 
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Figure  2  Point  charge  response  with  additive  white  noise,  SNR  =  5.11  dB 
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Figure  3  Input/output  relationships  for  the  mean  and  median  filter  of  length  three 
(L  =  3).  The  input  data  sequence  length  is  12  (N  =  12).  The  output  stays  equal  to 
zero  until  the  filter  is  completely  initialized. 
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Figure  4  Input/output  relationships  for  the  mean  and  median  filter  of  length  three 
(L  =  3).  The  input  data  sequence  length  is  1 12  (N  =  112).  The  output  stays  equal 
to  zero  until  the  filter  is  completely  initialized 
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Figure  5  Critical  portion  of  normalized  point  charge  response,  before  the 
explosion  point/first  zero  crossing.  The  miss  distance  is  equal  to  eight  meters 
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Figure  6  DFT  of  the  point  charge  response 
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Figure  7  Point  charge  response  with  additive  colored  noise,  SNR  =  -0.76  dB 


-30  -25  -20  -15 

Distance  (meters)  to  Zero  Crossing 


10 


15 


Figure  7.1   Mean  filtered  point  charge  response  with  colored  noise    The  noisy 
point  charge  response  is  presented  in  Figure  7. 
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Figure  8  Point  charge  response  with  additive  white  noise,  SNR  =  8.37  dB 
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Figure  8. 1    Combination  mean/median  filtered  noisy  point  charge  response    The 
noisy  point  charge  response  is  presented  in  Figure  8    The  threshold  point  is 
represented  by  (*),  and  the  zero  crossing/fuzing  point  is  represented  by  (o). 
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Figure  9  Point  charge  response  with  additive  white  noise,  SNR  =  4.86  dB 
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Figure  9. 1   Combination  mean/median  filtered  noisy  point  charge  response    The 
noisy  point  charge  response  is  presented  in  Figure  9.  The  threshold  point  is 
represented  by  (*),  and  the  zero  crossing/fuzing  point  is  represented  by  (o). 
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Figure  10  Reliability  test  of  the  combination  mean/median  filter. 
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Figure  11   Robustness  of  the  mean/median  filter  combination    Noise  variance 
between  0.01  (SNR  =  8.2086  dB)  and  0.02  (SNR  =  5. 1983  dB)    The  mean  error 
is  represented  by  (*),  one  standard  deviation  from  the  mean  error  is  represented  by 
(+),  100  trials  per  SNR  realization 
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Figure  12  Band  limited  Gaussian  noise  (colored  noise) 
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Figure  12. 1   Frequency  spectrum  for  the  colored  noise  presented  in  Figure  1 2    Fc 
is  the  center  frequency  in  Hertz,  fs  is  the  sampling  frequency  in  Hertz,  Q  is  the 
quality  factor,  and  N  is  the  order  of  the  FIR  filter. 
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Figure  13  Band  limited  Gaussian  noise  (colored  noise). 
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Figure  13.1   Frequency  spectrum  of  the  colored  noise  presented  in  Figure  13    Fc 
is  the  center  frequency  in  Hertz,  fs  is  the  sampling  frequency  in  Hertz,  Q  is  the 
quality  factor,  and  N  is  the  order  of  the  FIR  filter 
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Figure  14  Band  limited  Gaussian  noise  (colored  noise) 
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Figure  14. 1   Frequency  spectrum  of  the  colored  noise  presented  in  Figure  14    Fc 
is  the  center  frequency  in  Hertz,  fs  is  the  sampling  frequency  in  Hertz,  Q  is  the 
quality  factor,  and  N  is  the  order  of  the  FIR  filter 


Fc=1e+04    fs  =  4e+04    Q  =  8  and 


N  =  48 


0 


0.1 


0.2  0.3 

Frequency  (Hz) 


0.4 


0.5 


2X 


REFERENCES 

1  Hull,  D.  H.,  Internal  Report,  U.  S.  Army  Research  Laboratory,  Aldelphi,  Maryland, 

September  1993. 

2.  Strum,  R.  D.,  and  Kirk,  D.  E  ,  Discrete  Systems  and  Digital  Signal  Processing,  Addison 
Wesley,  1989. 

3.  Rabiner,  L  R  ,  Sambur,  M.  R  ,  and  Schmidt,  C.  E  ,  "Applications  of  a  Nonlinear 
Smoothing  Algorithm  to  Speech  Processing,"  IEEE  Transactions  on  Acoustics,  Speech, 
and  Signal  Processing,  VOL  ASSP-23,  NO  6,  pp.  552-557,  December  1975. 

4.  Fassnact  F.,  Principles  of  Systems  Engineering,  ECE  4010  Class  Notes,  U.  S.  Naval 
Postgraduate  School,  Monterey,  California,  Spring  1993. 


29 


Appendix  MATLAB  code  implementation 

%  GNF.m     John  K.    Dewey 

%  Generates   Colored  Noise, formed  as,    Band  Limited  Gaussian  Noise 

%  Uses    function   finddbf.m  to  measure   Center  Frequency  and 

%  Q  values. 

%  Matlab  toolbox  functions  used:  FIRl.m,  FREQZ.m,  SIZE.m 

clear  h   Clears  all  previous  variables 

%  Input  Arguments 

Fc  =  input ('Enter  the  center  frequency  =  ' ) ;  GO  =  num2str (Fc) ; 

Q  =  input ('Enter  the  value  for  Q  =  ');  Gl  =  num2str(Q); 

fs  =  input ('Enter  the  sampling  frequecy  f s  =  ' ) ;  G2  =  num2str(fs); 

%  Translate  Q  to  cutoff  frequencies 
Fl  =  Fc  -  Fc/Q/2; 
Fu  =  Fc  +  Fc/Q/2; 

%  Input  can  be  set  to  accept  the  cut  off  frequencies  directly 

%   Fc  = 

%    Fl  = 

%    Fu  = 

%   Q  -  Fc/(Fu-Fl) ; 

%  Specs 

yspecs  =  [-100  -100  -3  -3  -100  -100]; 

xspecs  =  [0  Fl  Fl  Fu  Fu  fs/2]; 

%  Normalize  Specs 
xnspecs  =  xspecs  . /fs; 

%  Set  warning  to  eliminate  aliasing 

if  fs/2  <  Fu  +  Fc/Q/2 

Warn  =  ['  ';'Warning  Aliasing  Warning']; 

disp(Warn) 

break 

end 

%  Find  the  order  of  the  filter 
N  =  28;  %  The  first  guess 
Qz  =  1;  %  The  minimum  Q 

while  Qz   <  Q 

F   =    [2*Fl/fs    2*Fu/fs]; 

[B]    =    firl(N,F); 

[H]  =  freqz (B, 1,2048) ; 
magz  =  20*logl0 (abs (H) ) ; 
xax  =  0: (fs/2) / (length (magz)-l) :fs/2; 

[WCz,Qz]  =  finddbf (magz, xax) ; 
N  =  N+2; 
end 
N  =  N-2;  G3  =  num2str(N); 
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%  Generate  White  Gaussian  Noise 

t  =  1:2048; 

wn  =  randn (size  (t ) ) ; 

%  Color  White  Noise 
Hf  =  filter (B, l,wn) ; 
S  =  spectrum  (Hf, 2  048) ; 
[n,m]  =  size  (S)  ; 
y  =  20*loglO(S(2:n)) ; 
x  =  (0:fs/2/(n-2) :fs/2) ; 
xn  =  x  . /fs; 

%  Display  the  Frequency  Spectrum  of  the  Colored  Noise 

plot (xnspecs, yspecs, xn, y) 

title (' Spectrum  of  the  Colored  Noise') 

gtext(['Fc  =  ',G0,'    fs  =  ',G2,'    Q  =  ',G1,'   and  N  =  ',G3]) 

xlabel ('Frequency  (Hz) '), ylabel ('Magnitude  (dB)'),grid 

figure 

%  Display  the  Colored  Noise 
plot (1 :length (Hf ) , Hf ) 
title ('  Colored  Noise') 
xlabel  ('  Index'  ) 
ylabel  ('Magnitude'  ) 
grid 

%  Save  the  Colored  Noise  for  future  use 
%save  en  Hf 
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%  FINDDBF.m,  Fall  1993,  John  K.  Dewey 

%  A  Function  that  measures  the  Center  Frequency  and  Q  values 

%  [wc,Q]  =  finddbf (mag,  f) 

%  input  magnitude  and  frequency 

function [WC,Q]  =  finddbf (mag, f) ; 

%  Find  the  Center  Frequency  Point 
[pk,  Ic]  =  max (mag); 

%  Find  the  -3dB  cut  off  frequencies 
A  =  abs (mag  -  pk  +  3); 
[xl,ql]  =  min(A(l:Ic)  )  ; 
[xll,qr]  =  min (A (Ic: length (A) ) ) ; 

%  Find  the  Center  Frequency 
WC  =  f (Ic); 

%  Compute  Q 

Q  =  f (lc)/((f (qr+lc)-f (ql))); 
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